Convexity of the self-energy functional in the variational cluster approximation 
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In the variational cluster approximation (VGA) (or variational cluster perturbation theory), 
widely used to study the Hubbard model, a fundamental problem that renders variational solutions 
difficult in practice is its known lack of convexity at stationary points, i.e. the physical solutions can 
be saddle points rather than extrema of the self-energy functional. Here we suggest two different 
approaches to construct a convex functional f2[E]. In the first approach, one can show analytically 
that in the approximation where the irreducible particle-hole vertex depends only on center of mass 
coordinates, the functional is convex away from phase transitions in the corresponding channel. 
Numerical tests on a tractable version of that functional show that convexity can be a nuisance 
when looking for instabilities both in the pairing and particle-hole channels. Therefore, an alterna- 
tive phenomenological functional is proposed. Convexity is explicitly enforced only with respect to 
a restricted set of variables, such as the cluster chemical potential that is known to be otherwise 
problematic. Numerical tests show that our functional is convex at the physical solutions of VGA 
and allows second-order phase transitions in the pairing channel as well. This opens the way to the 
use of more efficient algorithms to find solutions of the VGA equations. 

PACS numbers: 71.10.Fd, 71.15.Dx, 71.10.-w, 71.27.-fa 



I. INTRODUCTION 

Effective functionals, such as the Landau-Ginzburg 
free energy functional in the vicinity of phase transitions, 
have long been used to study classical and quantum sys- 
tems. Typically, the effective functional F\h\ is obtained 
in terms of a relevant variable h, such as the order pa- 
rameter in the Landau-Ginzburg theory, or the electron 
density in the density functional theory^. The optimal 
value of h is then obtained from the requirement that the 
functional be stationary at the solution, SF/Sh = 0. 

In the formalism of Luttinger WardS and Baym 
KadanofS, one of the best known functional approaches 
for correlated electrons, a functional fl[G] is stationary 
and equal to the grand potential at the exact physical 
value of G. In such a scheme, the functional depen- 
dence fl[G] is not known exactly. It can be approximated 
perturbatively by summing a subset of the infinite se- 
ries of skeleton diagrams that define the Luttinger- Ward 
functional. To obtain non-pcrturbative results on the 
other hand, one of the most effective methods is the 
Dynamical Mean-Field Theory (DMFT)i. Chitra and 
Kotliar— have shown how this theory can be obtained 
by modifying the Kadanoff-Baym functional and mak- 
ing the local approximation on the stationarity condi- 
tion. More recently, a general scheme for generating a 
wide class of non-pcrturbative approximations for the 
Hubbard model^ from a functional has been proposed 
by PotthofS. In this method, known as the self-cnergy- 
functional approach (SFA), a new functional of the 
self-energy S is constructed, which is stationary at the 
physical solution. The functional itself is unknown ex- 
plicitly, but Potthoff suggested a particular way of calcu- 
lating a variational solution Sft/SJ^ ~ with the help of a 
reference system, typically a cluster of finite size, which 
can be solved exactly. This particular implementation 



of the self-energy functional with no bath (contrary to 
DMFT) goes under the name of the variational clus- 
ter approximation (VGA). DMFT and generalizations 
thereof, (such as Gellular Dynamical Mean-Field The- 
ory (GDMFT)) can be obtained as various special cases 
of SFAfii^ corresponding to different choices of reference 
systems and/or approximations of the stationarity con- 
dition. (The functionals of Ghitra and Kotliar— and of 
Potthoffii^ are in fact identical, as shown in Appendix 

EH). 

One desirable feature of functional approaches is the 
variational principle that guarantees that an approximate 
grand potential is an upper bound to the true grand po- 
tential. Such a variational principle is missing for the 
stationary solutions of both the Baym-Kadanoff func- 
tional^ and Potthoff's self-energy functional^ since sta- 
tionary solutions are known to be saddle points rather 
than extrema. 

A question thus remains open up to present, whether 
or not it is possible to construct a functional, be it a 
functional of the Green function G or of the self-energy 
S, such that its stationary solutions would always be ex- 
trema (say, minima) of the functional. If so, this would 
mean that the functional is convex at the physical solu- 
tion, which is a first necessary step on the way to prove 
the variational principle. For an infinite-coordination 
Bethe lattice, this question was answered positively by 
Kotliar—, who proved that a functional could be con- 
structed, such that its extrema occurred at the physical 
local Green function of the Hubbard model. However an 
attempt by Chitra and Kotliar to find its analogue for a 
finite-dimensional lattice was only partially successful^. 

Another motivation to find convex functionals is a 
practical one. In VGA, the functional is definitely not 
convex for example when the intra-cluster chemical po- 
tential is varied. The physical solutions are then sad- 



die points. However, most efficient numerical algorithms 
such as e.g. the conjugate gradient method, have been 
designed to find extrema of a functional rather than sad- 
dle points. Although one may attempt to use such al- 
gorithms to also find saddle points (by minimizing the 
magnitude square of the gradients of the functionaU^), 
the unphysical solutions may occur. 

In this paper we reexamine the above problem and 
show that a convex functional can be found. In particu- 
lar, we construct a new functional of the self-energy r2[S], 
such that its stationary solutions are minima when the 
irreducible particle-hole vertex depends only on center 
of mass coordinates and the system is away from phase 
transitions in the corresponding channel. Going beyond 
this approximation involves complicated integrals that 
cannot be treated analytically, however numerical tests 
on a tractable version of the functional suggest that it is 
indeed always convex at the physical solution. Moreover, 
we show that such a construction is not unique and that 
several functionals can be constructed that differ in the 
higher-order terms of the expansion with respect to the 
sclf-encrgy. 

Despite the convexity of the proposed functional, its 
implementation requires additional approximations and 
it turns out to be inadequate to detect a second-order 
phase transition for both Cooper pairing and antifer- 
romagnetic instabilities. There, the new functional ap- 
pears to always be convex at the paramagnetic solution, 
whereas in the case of second-order phase transitions, the 
paramagnetic solution is rightly expected to be a saddle 
point, with minima developing instead at a finite value of 
the symmetry-breaking order parameter. In other words, 
the sought convexity of the functional 'overdoes' its job, 
imposing too stringent conditions on the resulting phys- 
ical solution. 

Given that the VGA method has originally been de- 
veloped to study broken symmetry phases, the aforemen- 
tioned feature of the proposed new functional is especially 
undesirable. To cure this drawback, we propose a differ- 
ent, this time phcnomcnological, approach that ensures 
convexity at the physical solution and, moreover, respects 
the tendency of the system to develop an instability to- 
wards a broken symmetry phase. This finding opens up 
new perspectives in the use of the convexity property of 
the functional, in particular permitting to apply power- 
ful numerical techniques, such as the conjugate gradient 
method, which are only guaranteed to work if the solu- 
tion is known to be an extremum (and not a saddle point) 
of the functional. 

The paper is organized as follows. First, in Section HIl 
we review the variational cluster approximation (VGA) 
as proposed by Potthoff?, followed by a general discussion 
of the stability of the stationary solution and criteria for 
convexity in Section lTTTl We then proceed to derive a new 
functional of self-energy in Section IIVI with the proof of 
its convexity given in Appendix |X] and the recipe for in- 
corporating it into the VGA framework in Appendix |B] 
The second part of Section |TV] is devoted to numerical 
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FIG. 1: (Color online) Luttinger-Ward functional ^[G] con- 
structed as a sum of renormalised skeleton diagrams with ap- 
propriate combinatorial coefficients (not shown) . The Hartree 
approximation consists in considering the first two diagrams 
only of the series. 

tests of the proposed functional on a two-dimensional 
Hubbard model, which show that the functional is convex 
but that it fails to correctly describe second-order phase 
transitions into a broken symmetry phase. Triggered by 
this negative result, we propose in Section |V] an answer 
to the convexity problem while correcting the aforemen- 
tioned failure to describe second-order phase transitions. 
The adequacy of this new construction is corroborated 
by numerical tests. More technical aspects of the work 
are detailed in appendices. 

II. REVIEW OF THE CONVENTIONAL VGA 
SCHEME 

Consider, following Potthoffiiii, a general Hamiltonian 
H = Holt) + Hi(\J) with one-particle hopping parame- 
ters t and two-particle interaction parameters U : 

H = J2 + ^ 51 U^jklc\c]clCk. (1) 

ij ijkl 

Here j, fc, I refer to an orthonormal and complete set of 
one-particle basis states. The equilibrium thermodynam- 
ics and elementary one-particle excitations of the system 
for temperature T and chemical potential arc fully de- 
scribed by the one-particle Matsubara Green functioriii^ 
defined by the imaginary-time ordered product T 

G.(ri,Ti;r2,T2) = -(T[c,(ri, ri)ct (rs, rs)]), (2) 

where the symbol (...) denotes thermal average 
and Gct(1,2) = Gcr(ri, ti; r2, T2) can be seen as 
(ri, Ti|GCT|r2, ^2) or G12, a matrix of two indices, each 
of which stands for both space and imaginary time. 

We start by defining the Luttinger-Ward functional 
<f>[G], constructed formally as a sum of all closed, irre- 
ducible skeleton diagrams involving fully renormalized 
("dressed") Green functions^, as illustrated in Fig.[T] We 
note that this functional is universal, that is, it docs not 
depend on the kinetic part of the Hamiltonian _ffo(t), but 
only on the vertices U and the Green function itself. 

The important property of the Luttinger-Ward func- 
tional is that its functional derivative gives the self-energy 
of the system: 

-M- = S.(2,1)==E.[G]. (3) 
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This equation serves as a definition of the self-energy as 
a functional of G. We can now introduce the grand po- 
tential of the system as the Baym-Kadanoff functional^ 
defined in terms of the fully dressed Green function as 
follows: 

^bk[G] = $[G] - Tr ((Go 1 - G'^) G) + Trln(G). (4) 

Inverting Eq. ^ (locally^) to obtain the Green func- 
tion as a functional of S, we can now express the 
Luttinger-Ward functional as a functional of self-energy 
<i>[G[S]]]. Substituting this into the expression for the 
Baym-Kadanoff functional ^ and using Dyson's equa- 



exactly (for example by means of exact diagonalization). 
The one-body part of the cluster Hamiltonian can con- 
tain a different hopping matrix, along with site energies, 
chemical potential and Weiss fields, all of which are used 
as variational parameters. Since one can write the grand 
potential for the cluster problem as 



= F'p'] - Trln [G'^^ - S') 



(10) 



the universality of the functional allows one to find 
its exact value for the solution of the cluster problem: 



Fp'] = = n' + Trln (Gq^^ - S' 



(11) 



tion G 



Gn 



E, Potthoff proposed the following 



functional of the self-energy 



Using this last result, the functional 17 [S] in Eq. (O can 
be evaluated exactly when S — > S' 



Op] = $ [G[q] - Tr (SG) - Trln(Go-i - S). (5) m'] = {^'[^'] + TrlnCG^^ - S')) - Trln(Go 

Recognizing that the first two terms on the right hand 
side of the Potthoff functional represent the Legendre 
transform of the Luttinger-Ward functional 



J^[S] = $ [G[S]] - Tr (EG) ^ 



(6) 



one can now rewrite the expression for the grand poten- 
tial as follows: 

0[E] = F[E] - Tr In {G^^ - S) . (7) 

it is easy to show that the fol- 



Using Eqs. ^ and % 
lowing equation holds 



E'). 
(12) 

Since they are the solutions of the cluster problems, 
the self-energies S' can be varied through the one-body 
terms of the clusters (which we define to also contain 
Weiss fields for various order parameters). These are 
collectively represented by the matrix t'^y Clearly, the 
self-energies obtained in this way will span only a small 
subspace of an infinite-dimensional space of all possible 
variations, namely only those that can be represented as 
the physical self-energies of a cluster E'(t') parametrized 
by the matrix t'^^ . The corresponding stationary solution 
is obtained by searching for values of such that 



Sl^a {1,2) 



= -G, (2,1) = -G2i[E] 



(8) 



an 



JE 



dE 



(13) 



It can be viewed as the definition of the Green function in 
terms of the self-energy E. Using this, we immediately 
arrive at the conclusion that the variational derivative 
of the grand potential in Eq. ([7]) vanishes at the true 
physical solution: 



50[E] 



ST, 



-G+ (Go 







(9) 



The set of equations (|12p - (fT3|) forms the essence of the 
VGA quantum cluster method. 



III. STABILITY OF THE STATIONARY 
SOLUTION 

Let us consider fluctuations around the stationary so- 
lution of an arbitrary self-energy functional: 



sol 



by virtue of the Dyson equation. At the physical solu- 
tion, 0[E] — il[G] is equal to the true grand potential. 
Therefore solving the problem amounts to finding such a 
function E that satisfies the above stationarity condition. 

Since the Luttinger-Ward functional $[G] is a univer- 
sal functional of the interaction U and of the stationary 
value of G, its Legendre transform F[E] is a universal 
functional of U and E. In other words the value of this 
functional should not depend on one-body operators such 
as the hopping matrix t. This last fact is the crucial point 
that allows one to define VGA. 

One proceeds as follows. In the case of the Hubbard 
model where the interaction is local, one can modify the 
hopping matrix elements to subdivide the infinite cluster 
into disjoint identical clusters. One can then compute the 
grand potential fl' and the self-energy E' of that problem 



Sn = 0[E + ST.] - n[T 



(14) 



s^n 



(l',l)5E„K2',2) 



JE,,(2',2) 



The spin indices have been written explicitly and must 
be summed over. From now on, we will not write them 
explicitly to have a lighter notation. 

For a stationary point to be numerically stable, it 
must be a minimum. A maximum will also do since it 
suffices to change the sign. Gorrespondingly, the func- 
tional derivative in Eq. p4p must be negative or positive- 
definite, respectively. It is easy to verify that for Pot- 
thoff 's self-energy functional, Eq. ([7|), the second func- 
tional derivative is given by 



s'-n 



STn'ST22' 



ril';22' + G112G2' 



(15) 
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where Tui -22' stands for the second functional derivative 
of the universal functional PI'S]: 



rii':2 



(5Sii'(5S]22' 



(16) 



and is thus a tensor of the fourth rank. 

It follows directly from the expression for the sec- 
ond functional derivative Eq. (fT5|) that the stationary 
solution is a saddle point. Following the arguments of 
Refi^, this is simplest to illustrate in the absence of two- 
body interactions in the Hamiltonian. when E = and 
= 0. Then the first term in Eq. (fT5|l vanishes while 
in Matsubara-Fourier space the last term leads to 



6^n 



= G {k[) G [k'^) 5 {k[ - fca) 5 [k'^ - fci) . (17) 



In the sector of zero total momentum and total energy 
{k'l = —k'2) , this quantity decouples in 2 x 2 blocks with 
zero diagonal elements and equal off-diagonal compo- 
nents equal to (w^ -I- ^^J"^- The eigenvalues are thus 
both positive and negative, which corresponds to a sad- 
dle point. The effect of interactions, Tiii-^22', cannot cure 
this problem for all wave vectors and frequencies. 

In the case of VGA, the relevant question concerns sta- 
bility with respect to variations in the cluster parameters. 
It has been pointed out empirically in the original VGA 
proposal^, that variations with respect to site energies or 
chemical potentials lead to saddle points. This is illus- 
trated in Fig. [2] where the grand potential obtained from 
a 2 X 2 cluster solution is plotted as a function of the 
intra-cluster nearest-neighbor hopping t' (Fig. [5^) and 
of the cluster chemical potential ^' (Fig. Wp)- Glearly, 
the stationary solution is a minimum in the former case 
and a maximum in the latter. Yet, varying the cluster 
chemical potential jj! has been showni^ essential in the 
VGA scheme for obtaining thermodynamic consistency, 
which requires that the number {n) of electrons in the 
system has the same value when calculated from the two 
independent relations: 



(n) 



9/i 



1 



Auj f{u) Im [TrG(w + iO+)] , 



(18) 

where /(w) = [exp(tj/T) -f 1] ^ is the Fermi function. 
Therefore, it is preferable to let ^' vary and somehow 
deal with the fact that the grand potential is known to 
be non-convex in this case. 

The fact that the solution is a saddle point rather than 
an extremum has consequences for practical implemen- 
tations, since there seem to exist no robust numerical 
search algorithms for a saddle point, whereas many such 
algorithms have been developed for extremaM searches. 
It would therefore be desirable from this point of view, 
as well as for the reasons outlined in Section |T1 to find 
such a functional 17 [E] whose stationary solutions would 
be guaranteed to be cxtrcma. In the next Section, we 
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FIG. 2: Functional Q, for the 2D Hubbard model {U /t = 8) 
calculated from the exact-diagonalization solution of a 2x2 
cluster using Eq. ( I12|) . plotted as a function of (a) cluster 
nearest-neighbour hopping parameter t' and (b) cluster chem- 
ical potential fj,' . The variational solution of the problem in 
the space of cluster parameters {t' = 1,^' = 4) is a saddle 
point instead of an extremum. 



shall demonstrate that such a functional can indeed be 
formally constructed, and its convexity can be proved 
rigorously given certain approximations. 



IV. CONVEX FUNCTIONALS n[S] 

In this section, we first derive a general expression for a 
convex functional. The drawback of this approach is that 
in practical implementations, convexity is preserved for 
order parameters as well. This makes it an undesirable 
feature in the presence of phase transitions. In the next 
section, we will restrict ourselves to preserving convexity 
as a function of the cluster chemical potential, which in 
practice is the main problem to be solved. 



A. Analytical results for convex functionals 

The main idea of the current approach is to add a term 
Af2[E] to the grand potential that is quadratic in E and 
that docs not alter the stationary solution of the original 
potential 12 [E], but renders it a positive-definite func- 
tional of the self-energy. We first introduce an auxiliary 
functional f[E] that is defined by 



f[E] 



f(l,2)E 

<SE (2,1) 



'5E(2,1) 



(19) 



S)"'(l,2) 



and that, by virtue of the stationarity condition Eq.Q, 
vanishes identically at the stationary solution of 12 [E]. 

It turns out that the simplest term A17[E] that can be 
added to the grand potential is of the form 



An, = -iTr([(Goi-E)f]' 



4- 



(20) 



5 



Indeed. Aril and its functional derivative both vanish at 
the stationary solutions of the grand potential, as follows 
from the definition of the functional f[E] and, therefore, 
the new functional 

r^ip] = rip] + Af7i[s] (21) 

yields the same stationary solution as the original one. 
The proof of the convexity of this new functional in the 
special case where the irreducible particle-hole vertex is 
local is given in the Appendix [X] 

We note in passing that our choice of Ar2[S] is not 
unique. To second order in the quantity f[S], the correc- 
tion Afli given by Eq. (|20p is the same as, for example, 
the following one 

AO2 = Trln (1 - {Ga^ - S)f) + Tr ((G^^ - E)f) (22) 
-Trln((G--I^)§)+Tr(l + (G--E)g), 

meaning that it leaves the stationary point of the original 
functional unchanged and leads to exactly the same ex- 
pression for the second functional derivative as Eq. (pOjl . 
Note also that despite being similar in spirit to the work 
of Chitra and Kotliar— , the above derivation is signifi- 
cantly different, as demonstrated in Appendix [Cl where 
the connection between the Chitra-Kotliar result and the 
Potthoff functional is exposed. 

While implementing Eqs. (^01 in practice, a prob- 
lem arises since Eq. for the auxiliary functional f[I]] 
contains a functional derivative 5i^[S]]/(5E whose value is 
not readily available in practical calculations. Therefore, 
an approximation must be made in order to implement 
the new functional into the VGA scheme. One such el- 
egant approximation is given in Appendix [B] and the 
conclusions of the numerical verification of its convexity 
is the subject of the following subsection. 

B. Drawback of this type of functionals 

While numerical tests show convincingly (see Ap- 
pendix |B] for details) that the functionals proposed in 
Eqs. ([20]) (jB5|) are indeed convex with respect to all clus- 
ter variational parameters, including the cluster chemical 
potential /i' , there is a major drawback. The numerical 
tests show one undesirable property, namely its failure to 
correctly predict second-order phase transitions, at least 
for the approximation to f[E] proposed in Appendix [Bl 

Consider Fig. [3^ that shows the dependence of il and 
of the the new functionals given by Eqs. (|B4p . on 
the magnitude of the antiferromagnetic Weiss field M. 
Unlike the original functional, which has a minimum at 
a finite value of M, both new functionals have a single 
minimum at i\/ = 0, thereby favouring the paramag- 
netic (PM) solution. The behavior as a function of the 
superconducting d-wave Weiss field, D, follows the same 
pattern, as shown in Fig.[3)D. Only the original functional 
n exhibits a minimum at non-zero value of D, whereas 
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FIG. 3: (Color online) Dependence of the original Potthoff 's 
functional Q (open circles), functional fii (solid red line) and 
Q2 (broken blue line) on the symmetry-breaking Weiss fields 
conjugate to (a) antiferromagnetic staggered magnetization, 
(b) d-wave superconductivity. Results are for the Hubbard 
model {U/t = 4). Arrows denote positions of minima of fl 
and maxima of fii. 

both functionals proposed in this work have no other min- 
ima except at D = Q. Clearly, such behavior of the new 
functionals is unphysical since the existence of antiferro- 
magnctism in the Hubbard model, for example around 
[/ = 4 at half-filling, is solidly established. 

As discussed in Appendix [XI we would have expected 
that at least in the particle-hole channel, instabilities 
(such as antiferromagnetism) could have been detected 
by the present functional. However, we had to approxi- 
mate (5i^[S]]/(5I] to do the calculation. In some sense, the 
requirement of convexity here 'overdoes' its job by ren- 
dering the paramagnetic solution always a minimum and 
hence ignoring a possible instability towards a broken- 
symmetry phase! 



V. CURING THE PROBLEM OF CONVEXITY 
ARISING FROM A RESTRICTED SET OF 
VARIABLES 

In practical implementations of VCA, the cluster 
chemical potential leads to a first-order saddle point even 
in paramagnetic states. There exist several numerical 
techniques that can deal with this type of problemsi^. 
In the context of VCA we show a practical way to trans- 
form the saddle point coming from variation of the cluster 
chemical potential /x' into a minimum while preserving 
the 'right' of the system to develop a second-order phase 
transition. It seems to work well even in cases where the 
analytically-obtained approximations given in Eqs. (|B4[) . 
(lB5l) fail. 

Let us recall that the general form of the convex cor- 
rection developed in Section [IV] was 

Typically, we know from experience about the existence 
of a certain variational parameter h (here, the cluster 
chemical potential /i') such that the fl lacks convexity 
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FIG. 4: (Color online) Dependence on the cluster chemical 
potential /x' of Potthoff's Q, (black dots), and the functional 
il\ in Eq. (|25|l for two choices of the coefficient A: the constant 
A = 15 (solid red line) and the variable A = 2Ao(/i') (dashed 
blue line). Results for the half-filled Hubbard model with (a) 
U/t = 4 and (b) U/t = 8 are shown. 



at the stationary solution h = Hq. 
correction in Eq. as follows: 



Let us modify the 



An 



'sn as 



/on 



(24) 



where we have effectively substituted the unknown func- 
tional derivative SQ/ST, by a much simpler derivative 
with respect to the variational parameter /i, which can 
be easily calculated numerically. 

We thus postulate the following functional 



(25) 



where A is some empirical coefficient that should be cho- 
sen such that the resulting potential Q\ is a convex func- 
tion of h. By construction, the additional term vanishes 
at the stationary solution h = where dil/dh = 0, and 
the value of the new potential coincides with the old one 
ri(/io)- For the second derivative at ho we have: 



h—hn 



'dh? 



1 + A 



'dh? 



(26) 



If the original hmctional is already convex, d^il/dh^ > 0, 
the new functional will be convex too. If however Q has 
a maximum and not a minimum at ho, the coefficient A 
has to satisfy the following inequality to ensure that the 
new functional is convex: 



A > Xo{h) = 



\dh^ ) 



(27) 



In particular, note that choosing A — 2Xo{h) will result 
in the new functional being convex everywhere and hav- 
ing the same absolute value of curvature as the original 
one: d^il\/dh^ = —d^fl/dh^ . Chosen in such a way, 
the behavior of the newly constructed functional fix as 
a function of the cluster chemical potential h = fi' is 
plotted in Fig. g) 

As we see, fix can clearly be made convex with a suit- 
able choice of the coefficient A. The lowest allowed value 



Ao that yields a convex functional is itself a function of 
parameters of the model, such as Hubbard U, as is clear 
from comparison of Fig. 3^) and b). In particular, we 
find Ao = 5.8 for U = 4t and Ao = 9.9 for U = 8t. 

Note that by taking derivatives with respect to clus- 
ter parameters in Eq. (|24p instead of the more general 
Eq. (j23p . we have substituted a strong requirement of 
convexity with respect to all variations in S with a much 
weaker one, which only requires that the functional be 
convex with respect to a particular parameter (or set of 
parameters) h. Let us then see what effect this has on 
the dependence of J^a on a symmetry-breaking parame- 
ter, such as the superconducting Weiss field D. The first 
derivative with respect to D reads: 



dnx _ dn d^n on 

ao ^ ao ^ dhdD~dh 



(28) 



All the stationary points of (where derivatives with 
respect to both /i' and D vanish) are also stationary 
points of the new functional fix- That functional can 
also have additional stationary points that are not sta- 
tionary points of O, but with the procedure suggested 
below to choose A, we found this not to be an issue. 

In numerical calculations, we searched for minima of 
VLx and found that, indeed, one recovers as minima the 
same solutions (/Uq, Do) that were saddle points of fi. To 
illustrate how one can choose A in practice, we consider 
several cases. First, in Fig. [5^ ^' is fixed at the known 
superconducting solution /ig and A is taken either as con- 
stant or defined by the procedure in Eq. (^7)) that guar- 
antees convexity. All curves have their minimum at the 
same value D = Do as expected. To gain more insight 
into the convergence process, we can also check whether 
the solution for D is stable in cases where the value of 
the cluster chemical potential jj! is slightly off the true 
solution n'o- Such a situation is illustrated in Fig. [5)3 
where the chemical potential ^' is fixed at the PM value 
/ipjj-j = 2.45 instead of the true SC solution ^'o = 2.28. 
As expected from Eq. (PS)) . both the old and the new 
functionals coincide at the paramagnetic solution D = 
since by construction, dfl/dfi' = there. The situation 
is different however away from D — 0. Fixing A = const 
(solid line in Fig. [SJj) yields a minimum at an incorrect 
value of D = 0.84 instead of Do = 0.56. By contrast, a 
variable coefficient A = 2Ao(/i'), where Ao is determined 
at each point from Eq. (|27p . gives a minimum of Qx{D) 
(dashed line in Fig. [SJd) that is already very close to the 
true solution. Varying fi' as well eventually leads to the 
correct solution in all cases, as mentioned above. An em- 
pirical case-by-case analysis has shown that the choice 
A = 2Xo{fi') gives consistently reliable results and is pre- 
ferred over fixing A to a constant value. Besides, such 
a choice avoids the ambiguity in the value of A and en- 
sures that the functional is always a convex one, 
according to the inequality in Eq. ((27| . 

The knowledge of the derivatives dil/dh and d^ft/dh^ 
is required in order to implement the convex correc- 
tion ([^5]) with A ~ 2Xo{h). Numerically, this requires 
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FIG. 5: (Color online) Dependence on the symmetry-breaking 
superconducting Weiss field D of Potthoff 's grand potential 
(black dots), and of the functional in Eq. p5|l constructed for 
two choices of the coefficient A: the constant A = 15 (solid red 
curve) and the variable A = 2Ao(/i') (dashed blue curve). The 
calculations are for the electron-doped Hubbard model (2.2% 
doping) with U = At and the next-nearest neighbour hopping 
t' — 0.3t. The cluster chemical potential fi' was fixed at the 
value corresponding to (a) the true SC solution, D — 0.56, 
and (b) the PM solution, D = 0. 

the knowledge of the functional at three points. For 
this reason, the computational cost of evaluating £7 at a 
given point in parameter space is three times that of the 
original functional Q. However, the minima of the pro- 
posed functional may now be searched more efficiently 
using powerful numerical methods designed for extrema 
search, such as the conjugate gradient methodic. 



is varied. Saddle points are notoriously more difficult 
to find numerically than extrema. We have thus con- 
structed a new functional that we proved to be a convex 
functional of the self-energy at the stationary solution, 
at least in the case where the irreducible particle-hole 
vertex depends only on the centre of mass coordinates 
(Appendix E]) . It turns out, however, that implementing 
the proposed functional in practice is far from simple, 
since it involves, as can be seen in Eq. (j20p . an unknown 
functional derivative 6F/5Y1 of the Legendre transform 
F[S\ of the Luttinger-Ward functional. An approxima- 
tion therefore had to be made to express this functional 
derivative in terms of cluster-defined quantities, as de- 
tailed in Appendix [B) The corresponding numerical re- 
sults in Sec. IIV Bl show that the new functional is in- 
deed convex even in the general non-perturbative case. 
However attractive this may be, the proposed functional 
changes qualitatively the behavior of Q with respect to 
symmetry-breaking order parameters, such as magnetism 
or superconductivity, rendering instead the paramagnetic 
solution stable. 

To cure this problem we proposed a functional Eq. (j25p 
and a (non-unique) recipe to find the coefficient A that 
guaranties convexity with respect to cluster chemical po- 
tential /i' only. This approach removes the saddle point 
normally associated with /i', has the same physical solu- 
tions as the original problem and leaves unchanged the 
minimum or maximum character of the functional with 
respect to variations of symmetry-breaking order param- 
eters, such as magnetism or superconductivity. 

Despite the non-uniqueness of the proposed convex 
functional and the admittedly phenomenological basis of 
its derivation, we argue that this is an important result. 
In particular, our findings open a way to the use of pow- 
erful numerical algorithms for solving for minima, such 
as e.g. conjugate gradients, that only work provided the 
functional is convex at the physical solution. We consider 
the use of such efficient algorithms as highly desirable for 
the VGA approach. The question of the existence of a 
functional that would give a bound for the true grand 
potential and would thereby implement the variational 
principle, is left open. 

We thank G. Kotliar for useful comments on the 
manuscript. Gomputations were performed on the Elix 
and the Ms RQGHP clusters. The present work was 
supported by FQRNT (Quebec), NSERC (Ganada), CFI 
(Ganada), CIAR, the Tier I Ganada Research Chair Pro- 
gram (A.-M.S.T.). 



VI. CONCLUSIONS 

It is known that physical solutions obtained with VGA 
are in general saddle points rather than extrema of the 
functional fl[Y,]. This is so in particular in the impor- 
tant practical case where the cluster chemical potential 



APPENDIX A: CONVEXITY OF THE 
ANALYTICAL FUNCTIONAL 

We aim to prove the convexity of the new proposed 
functional fli[Ti\ given by either of Eqs. (PO)) {i ~ 1) or 
dm (i = 2). At the true solution where Sfl^/ST, = f[E] = 
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FIG. 6: Bethe-Salpeter equation for the particle-hole vertex 
n in the language of Feynman diagrams. 



definitions 



5h2 

6n 



34 



S^-F ^ SGi2 

r21;34 + G13G42. 



where r2i;34 has the symmetry r2i;34 
tain, 



(A2) 
r34;2i, and ob- 



0, we have that 



22' 



f=0 



-Tr 



22' 



(Al) 



with G ^ ~ {Gq ^ — 5]) and matrix multiplication implied 
for the indices not explicitly written. We then use the 



l':22' 



^i^a ^ ^ S)33' • r43':ll' • (Gq — 2)44/ • r34';22 



(A3) 



f=0 



with the summation implied over the repeated indices. 
The last term may be rewritten in the following short- 
hand matrix notation: 



(Go 1 - E)33' • r 



43';11' 



(G, 



-1 




TT[G-'Tl,,G-'rl,,] 



(A4) 



r 



where the transpose applies only to the indices involved 
in the matrix mutiplication. 

Irreducible vertices are usually easier to approximate. 
So it is useful to work with the second functional deriva- 
tive of the Luttinger-Ward functional <&[G] with respect 
to the single-particle Green function: 



rph 



S^^G] 

SGyiSG; 



6Y. 



11' 



2'2 



<5G2' 



2', 2 



(A5) 



where we have used the property Eq. ([3]) of the Luttinger- 
Ward functional. The superscript in the notation 

indicates that this symmetric matrix (t 

is the irreducible particle-hole vertex, i.e 



12;34 — 34;12 



nr2' 



34 



(^13(5; 



I3O24 



G22'r 



ph 



22'-L 2'3';34 



G3'l, 



(A6) 



where 11 is the full (reducible) particle-hole vertex. In the 
language of Feynman diagrams, the latter can be cast as 
shown in Fig. [5) 

It is required on physical grounds that the reducible 
particle-hole vertex 11 should be positive-definite at zero 



frequency for stability in the particle-hole channel and 
hence it follows that its inverse obeys the same prop- 
erty. If there is an instability that affects the particle-hole 
channel directly or indirectly, it should still be visible in 

Using the identity 5F/5Yj = — G and the definition of 
r given by Eq. (|16p . we obtain (with summation over 
repeated indices implied) that 



p pph 
-L 12;1'2'-1- y 



<5(-G)21 (^Ei'2' 



2';34 



dLi/2' 0G43 

= -(52,4<5i,3 = (A8) 
which again has the structure of matrix multiplication if 
the pairs of indices on either side of the semi-colon are 
flattened (combined as one). Using this result we obtain 



Tr(G 



-It-T 



T\yG 



-IpT 



a2'mGT^^,.GTll,.) = (A9) 

'5l,3<5l',3' = I- 



-^ph 



Using this last identity and Eq. 
find: 



H for 5^n.,/5Y?, we 



J 



i5Eii'(5S22' 



Tr(G^P2^.G^g,.) - ,5i,3(5i',3' - Gi4rP^.33,G4'i' 



(AlO) 



As has been pointed out above, the right-hand side of this from the requirement that the kernel of the particle-hole 
equation is positive-definite at zero frequency, following 
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vertex equation be positive. 
Let us now concentrate 



Tr(GrP2,.Gr: 



33' 



on L22':33'[G] 

appearing on the left hand side of 
the last equation (|A10|) . We assume that the irreducible 
particle-hole vertex is local (like in the Random Phase 
Approximation), i.e. 



i 11/. 22' — Cii/022'i 12- 



(All) 



With this approximation, we can write for the left-hand 
side of the stability equation (jAlOp : 



(SSll(5S22 22;5'5'"^5'6'^ 6'6';33 



-^ph 



Taking the Fouricr-Matsubara transform, we find 



(5Eii(5S22 

(5Sii(5S22 



r^5/G6'5'G5'6'r^,''3 



j23 



/ ^G(k,iw„)G(k + Q 



-rpi^ iQ,iu;,) X (Q, «^.) rp'^ (Q,iw.) 

I 



(A12) 
(A13) 

(A14) 
(A15) 



where the dressed Lindhard function x(Qi i^y) (no ver- 
tex correction) denotes the value of the integral (defined 
with the minus sign) in the last equation. Using the 
spectral representation one can show that x(Q,*w,y) for 
a system at equilibrium is always positive for all Mat- 
subara frequencies. Also, T^^ {Q^ioj^) is real, as follows 

from its spectral properties, so (pph (Q ,iuj^)j is posi- 
tive. It is thus clear that the quantity L{Q, itu^) must be 



negative-definite, with possible exception of phase tran- 
sitions that, in this simple approximation, could appear 
in the particle-hole channel. Comparing this with the 
identity given by Eq. (|A10p : 



L[G] = 1 - GPP'^G, 



(A16) 



{Q,tLU^) ■ L (Q,iw^) = 1 + X (Q, iu;^) PP*^ (Q, iu;^) 



(A17) 



and remembering that its right hand side is positive- 
definite, we arrive at the conclusion that, at least in the 
above approximation for the particle-hole vertex, the ob- 
ject S'^rii/ST,'^ must also be positive-definite at zero fre- 
quency. This in turn means that our proposed grand 
potential defined by Eq. (|21|1 and either expres- 

sions ([20]) or ([22]) , is a convex functional of the self-energy 
within this approximation. 

It should be noted that being a convex functional of 
the self-energy is not necessarily the same as being a con- 
vex functional of a given cluster parameter h' . Indeed, 
differentiating again the first derivative dfl/dh' given by 



Eq. (|T3|) . one obtains 

d^fti _ rr s^n, dE(2)dS(i) r sn, d^Y.{i) 

'dhP ~ J J SJ:{l)6J:{2)~dh' dhJ~^J (5S(1) dh'^ ' 

12 1 

Note that in the above, we have assumed translational in- 
variancc, and hence self-energy can be written as a func- 
tion of one momentum and frequency variable only. Since 
the functional derivative in the first term of Eq. (jAlSp is 
almost certainly positive-definite at the stationary solu- 
tion, as suggested above, the first term is guaranteed to 
be positive too. The situation with the second term is 
more complicated. Naively, it may seem that this term 
vanishes since, at the stationary solution, SQ/STi is zero 
by definition. This would be true if the variational space 
of the cluster parameter h' was sufficient to describe the 
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complete variational space of the functional ri[E(/i')]. 
Unfortunately, and this is the main approximation en- 
tering the VGA method, this may not always be true. In 
fact, one may only vary ST,(h') in a subspace that can 
be parametrized by the cluster parameter h' , which is a 
small part of the whole infinite-dimensional space of vari- 
ation (5S. Therefore, the solution obtained by requiring 
that the derivative dfl/dh' vanishes, is not necessarily the 
same as the true solution where 6fl/6Y^ = 0. Although 
we expect this contribution to be very small, it is difhcult 
to judge the convexity of the given grand potential with 
respect to some cluster parameter, even if the potential 
is known to be a convex functional of S. Since varying h' 
is the best one can achieve within the VGA approach we 
stress that numerical tests are always desirable to verify 
the convexity of a given functional ri[I](ft,')]. 



APPENDIX B: IMPLEMENTATION OF THE 
CONVEX FUNCTIONAL IN VCA 

In Section Hv] we proposed a correction AfJi 2 to the 
original grand potential such that it is rendered a 

convex functional of the self-energy. The issue that we 
address here is how to implement this correction in the 
framework of the VGA quantum cluster method. 

The key element entering the expressions pO)) and ([22]) 
for Afi is the auxiliary functional f[S] defined by Eq. (fT9|) 
as a functional derivative f[E] = SH/S'S. Unfortunately 
it is not possible to evaluate this functional derivative 
even numerically since we have no direct way of vary- 
ing the self-energy S (however we develop a variant of 
this idea further in Section |V| . The unknown clement 
in f[S] Eq. ^ is the functional derivative 6F[i:]/ST., 
where F[Y^] is a universal functional of the self-energy, as 
explained in section |TT1 This functional derivative can be 
evaluated at the ground state solution of the cluster by 
virtue of Eq. ©: 



<5f [S] _ 6F[Y.'] 
(5E 



(Bl) 



sol' 



where, as before, we denote the quantities belonging to 
the cluster by prime. Now, that both terms on the right 
hand side of Eq. p9p are known, we can write: 



fp] ~-G' 



1 



Gn 



-G' + G. 



(B2) 



We stress that this last equation is approximate. This 
is because the self-energy E' of the cluster solution, at 
which the derivative in (jBl[) is evaluated, is not, generally 
speaking, equal to the true lattice self-energy. Another 
way to see this is to note that, by construction, f[S] must 
vanish identically at the true stationary solution of the 
lattice, whereas it is clear that the right hand side of 
Eq. (jB2[) can only vanish if the two Green functions are 
equal to each other at the ground state of the cluster 
tiling. Rigorously speaking, this can only be the case in 
the limit of infinitely large cluster, where G' — > G. 




FIG. 7: (Color online) Comparisons of three different func- 
tionals for the 2D Hubbard model ([//fiat ~ 8): the original 
functional from Eq. (|12|) (black dots), and the two function- 
al proposed in this work, given respectively by Eq. ()B4|l 
(solid red curve) and Eq. (|B5P (dashed blue curve). The 
panels show (a) dependence of all three functionals on the 
nearest-neighbour hopping parameter t' of the cluster, (b) 
same for the dependence on the chemical potential n' of the 
cluster. Comparison between the two functionals proposed in 
this work as functions of the cluster chemical potential n' for 
(c) half-filled case, /iiat ~ 41, and (d) away from half-filling, 
Miat ~ 3t (corresponds to 0.4% electron doping). 



Using this approximation Eq. (|B2|) , we can now rewrite 
Eq. ^ for Afli as 

An.m = -^Tr [1 - (Go-i - E)G')]' , (B3) 

so that, by virtue of Eq. (fT^ . the new convex functional 
ill = + Afii finally becomes: 

r2i[E]«f7' - Trln [(G(7^ - E)G'] 

- iTr[l-(Go-i-E)G')]'. (B4) 



Similarly, for the correction Ari2 expressed by Eq. (|22p . 
we obtain an alternative expression for the new functional 

n2p] + Tr [1 - (G(7^ - E)G'] . (B5) 



The proposed functionals Eas. (|B4p and (|B5p have been 
implemented into the VGA scheme using the Lanczos al- 
gorithm of exact diagonalization (ED) to solve the cluster 
problem. The Hubbard model on a square lattice with 
nearest neighbor hopping t has been studied, and a clus- 
ter of 2 X 2 sites was used in the ED scheme. While the 
hopping parameter t and the chemical potential /x of the 
lattice model remain fixed, we have the freedom of vary- 
ing the corresponding parameters of the cluster t' and 
fi' , whose optimal values should be found by solving the 
stationarity equation ([Te 
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Figure [7^ shows the comparison between the original 
grand potential Q (dotted line) and the two potentials 
derived in Section as functions of the variational clus- 
ter hopping parameter t' . All three functionals appear to 
have minimum at the same value of t' = t = I and the 
proposed new functionals are convex functions of t' near 
the solution, as is the original grand potential Q. 

The dependence of the grand potential on the cluster 
chemical potential fi' is plotted in Fig. [7)3 for the half- 
filled case. Functional Oi shown by the solid line, clearly 
is a convex function near the stationary solution fi^ = 
4:t, unlike the original grand potential fl (dotted line) 
that develops a maximum at this point. The details of 
the behavior of the new functionals near /j,q appear more 
clearly in the blow up shown in Fig. [71:. The functional 
^2 (dashed line) has vanishing second derivative within 
machine precision around the stationary solution. The 
same conclusion equally holds away from half-filling, as 
illustrated in Fig. [TJi. 

The above illustrations as well as many tests performed 
for different values of the parameters U, t and n all show 
that the proposed functional fii given by Eq. (|B4p devel- 
ops a minimum at the stationary solution with respect to 
both variational parameters ^' and t' . The functional fl2, 
is less useful for practical calculations since it exhibits a 
vanishing second derivative with respect to ^' sufficiently 
close to half-filling. 

Despite these encouraging results, note that the actual 
values of f2i and Q2 at the solution, while being nearly 
equal, deviate appreciably from the original grand po- 
tential f2. There would be no such deviation if we could 
have evaluated f[S] exactly instead of approximately as 
we were forced to do. 



where ^[G] is the Luttinger-Ward functional. Since 
J = when the Dyson equation is satisfied, both func- 
tionals are equal to the grand potential at the stationary 
solution. Chitra and Kotliar have shown^ that their new 
functional has a different stability criterion for its station- 
ary solution from that of the Baym-Kadanoff functional, 
however they have been unable to prove its convexity at 
the stationary point. In fact, they have shown explicitly 
that, in the Hartree approximation, their proposed func- 
tional ricK is unstable for repulsive interactions. This 
prohibits the use of the Chitra-Kotliar functional as a 
variational free energy and leaves open the question of 
convexity. 

Is there a way to relate the Chitra-Kotliar functional 
Eq. (jCip to the functionals of self-energy rii[S] and il2[S] 
proposed in Section IIVI ? From the previous two equa- 
tions, it is easy to show that (see Eq. (12) in Ref.— ): 



(5$ 

f^CK[G] = $[G]-Tr ( G— ) -Trln 



^" 5G 



(C3) 



Using the relation ([3]), S^/5G = S[G], and the Legendre 
transform of the Luttinger-Ward functional- as in Eq. ([6|) 



F[E] = $ [Gp]] - Tr (EG) , 



(C4) 



ficK in Eq. (|C3|) becomes 



APPENDIX C: CONNECTION TO THE 
CHITRA-KOTLIAR FUNCTIONAL 

The purpose of this Appendix is to identify the connec- 
tion that exists between Potthoff's original functional^, 
the new convex functional 51 [E] proposed in this work 
(Sec. IIVP , and an earlier attempt to construct a convex 
functional undertaken by Chitra and Kotliar in Refi^. 
According to the latter study, one can construct an im- 
proved version of the Baym-Kadanoff functional bk [G] 
given by Eq. ^ as follows: 

r!cK[G] = nBK[G] - Trln(l + ,/G) + Tr(JG), (CI) 

where J is the external source field coupled to the elec- 
tron's Green function 



11ck[G(E)] = F[Y] - Trln [Gq ^ - E] (C5) 



This last expression is nothing else but Potthoff's func- 
tional ri[E] as given in Eq. ([7]). In other words, we see 
that the Chitra-Kotliar functional of G, if expressed in 
terms of the self-energy E, is equivalent to Potthoff's 
ri[E]. Although both groups begin with the Baym- 
Kadanoff functional to propose their own functional, they 
use the Dyson equation in apparently different but in fact 
equivalent ways. 

As regards the new functionals of the self-energy intro- 
duced in this work (see Section [TV)) . their derivation, de- 
spite being similar in spirit to that of Chitra and Kotliar—, 
have required additional ingredients that are contained 
neither in the Potthoff's grand potential ri[E], nor in the 
Chitra-Kotliar functional. 
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